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ABSTRACT. We consider the design of an effective and reliable adaptive finite element 
method (AFEM) for the nonlinear Poisson-Boltzmann equation (PBE). We first examine 
the two-term regularization technique for the continuous problem recently proposed by 
Chen, Hoist, and Xu based on the removal of the singular electrostatic potential inside 
biomolecules; this technique made possible the development of the first complete so- 
lution and approximation theory for the Poisson-Boltzmann equation, the first provably 
convergent discretization, and also allowed for the development of a provably conver- 
gent AFEM. However, in practical implementation, this two-term regularization exhibits 
numerical instability. Therefore, we examine a variation of this regularization technique 
which can be shown to be less susceptible to such instability. We establish a priori es- 
timates and other basic results for the continuous regularized problem, as well as for 
Galerkin finite element approximations. We show that the new approach produces reg- 
ularized continuous and discrete problems with the same mathematical advantages of 
the original regularization. We then design an AFEM scheme for the new regularized 
problem, and show that the resulting AFEM scheme is accurate and reliable, by proving 
a contraction result for the error. This result, which is one of the first results of this type 
for nonlinear elliptic problems, is based on using continuous and discrete a priori L°° 
estimates to establish quasi-orthogonality. To provide a high-quality geometric model 
as input to the AFEM algorithm, we also describe a class of feature-preserving adaptive 
mesh generation algorithms designed specifically for constructing meshes of biomolec- 
ular structures, based on the intrinsic local structure tensor of the molecular surface. All 
of the algorithms described in the article are implemented in the Finite Element Toolkit 
(FETK), developed and maintained at UCSD. The stability advantages of the new reg- 
ularization scheme are demonstrated with FETK through comparisons with the original 
regularization approach for a model problem. The convergence and accuracy of the 
overall AFEM algorithm is also illustrated by numerical approximation of electrostatic 
solvation energy for an insulin protein. 
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1. INTRODUCTION 

The Poisson-Boltzmann Equation (PBE) has been widely used for modeling the elec- 
trostatic interactions of charged bodies in dielectric media, such as molecules, ions, and 
colloids, and thus is of importance in many areas of sciences and engineering, including 
biochemistry, biophysics, and medicine. The PBE provides a high fidelity mean-field 
description of the electrostatic interactions and ionic distribution of a solvated biomolec- 
ular system at the equilibrium state, and entails singularities of different orders at the 
position of the singular permanent charges and dielectric interface. The popularity of 
the PBE model is clearly evidenced by the success of software packages such as APBS, 
CHAPvMM, DelPhi, and UHBD. We summarize the mathematical PBE model in some 
detail in Section 2, referring to the classical texts [36, 48] for more physical discussions. 

While tremendous advances have been made in fast numerical solution of the PBE 
over the last twenty years (cf. [35, 25, 26] for surveys of some of this work), mathemat- 
ical results for the PBE (basic understanding of the solution theory of the PBE, as well 
as a basic understanding of approximation theory for PBE numerical methods) were fun- 
damentally unsatisfying, due to the following questions about the PBE and its numerical 
solution which remained open until 2007: 

1) Is the PBE well-posed for the dimensionless potential u! 

2) What function space does the solution u lie in? 

3) Can one derive a priori (energy and/or pointwise) estimates for the solution ul 

4) Is there an efficient (low-complexity) and reliable (provably convergent under uni- 
form mesh refinement) numerical method that produces an approximation Uh to 
the ul 

5) Is there a provably convergent adaptive method for the PBE? 

That these basic questions were open through 2007 is somewhat remarkable, given the 
popularity of this model. However, four key features of the PBE model, namely: (1) 
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the undetermined electrostatic potential at the boundary of a given system; (2) the singu- 
lar fixed charge distribution in biomolecules; (3) the discontinuous dielectric and Debye 
constants on the irregular dielectric interface (with a possible second interface represent- 
ing an ion exclusion layer); and (4) strong nonlinearity in the case of a strong potential 
or heavily charged molecules, place the PBE into a class of semilinear partial differential 
equations that are fundamentally difficult to analyze, and difficult to solve numerically. 
In fact, numerical evidence suggested that the most popular algorithms used for the PBE 
were actually non-convergent under mesh refinement, which would put the reliability of 
scientific results based on numerical solution of the PBE in doubt. 

To address this issue, in 2007 Chen, Hoist, and Xu [11] used a two-scale decomposi- 
tion as a mathematical technique to answer each of the open questions above about the 
PBE, building the first available solution theory and approximation theory for the PBE. 
(A basic existence and uniqueness result using variational arguments had appeared al- 
ready in [23].) A splitting-type treatment of the singular charges was not new, and is a 
very natural physical idea first sketched out in [21] and then also explored numerically 
in [64]. This method decomposes the PBE into a Poisson equation with singular charge 
and uniform dielectric that determines a singular function u s , and a regularized Poisson- 
Boltzmann equation (RPBE) that determines a smooth correction u, with the sum of the 
two giving the dimensionless potential: u = u s + u. This natural splitting technique 
was exploited in [11] to show that: the regularized PBE is well-posed, as is also the full 
PBE; the solution u can be split into a singular function u s (having a simple closed form 
expression) and a smooth remainder u which lies in a well-understood function space 
H 1+a with a ^ 0; the remainder function u is pointwise bounded almost everywhere; 
a standard finite element discretization that incorporates the singular function converges 
and does so at optimal rate in the limit of uniform mesh refinement; and finally, an im- 
plementable adaptive algorithm exists that can be proven mathematically to converge to 
the exact solution of the PBE. 

While this two-scale decomposition made a number of basic mathematical results pos- 
sible in [1 1], the resulting numerical algorithms (both based on uniform mesh refinement 
and adaptive mesh refinement) are subject to a hidden instability. This instability is in fact 
tied to the two-scale decomposition technique itself. (Example 2.1 in Section 2.2 gives 
a more complete description of this difficulty.) While this feature of the splitting has no 
impact on the mathematical or convergence results in [11], the practical impact is that 
algorithms based on this particular decomposition are not accurate enough to be compet- 
itive with other approaches. A slightly modified decomposition scheme was proposed 
by Chern et. al [14] (see also [20]) and was applied together with Cartesian grid-based 
interface methods to solve the PBE for simple structure geometries. The new splitting 
technique gives rise to a modified form of the regularized PBE with similar structure to 
the splitting scheme in [1 1], but appears to be more stable. 

This article is focused on using a similar decomposition variant to remove the insta- 
bility present in the formulation appearing in [11], as well as to improve the theoreti- 
cal results and algorithm components of the adaptive finite element algorithm described 
in [11]. In particular, we adopt a variation of the regularization splitting scheme similar 
to [20, 14, 1 1], involving a 3-term expansion rather than a 2-term expansion. We establish 
several basic mathematical results for the 3-term splitting, analogous to those established 
in [11] for the original 2-term splitting. This includes a priori L°° estimates, existence 
and uniqueness, and discrete estimates for solutions, and basic error estimates a general 
class of Galerkin methods. We then focus specifically on finite element methods, and 
design an adaptive finite element method (AFEM) for solving the resulting regularized 
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PBE. Due to recent progress in the convergence analysis of AFEM for linear and non- 
linear equations [10, 27], we also substantially improve the AFEM convergence result 
in [11] to one which guarantees contraction rather than just convergence. We present 
numerical examples showing the accuracy, efficiency, and stability of this new scheme. 
To provide a high-quality geometric model as input to the AFEM algorithm, we will 
also describe a class of feature-preserving adaptive mesh generation algorithms designed 
specifically for constructing meshes of biomolecular structures, based on the intrinsic 
local structure tensor of the molecular surface. 

While we focus on (adaptive) finite element methods in this article, the splitting frame- 
work we describe can be incorporated into finite difference, finite volume, spectral, 
wavelet, finite element, or boundary element methods for the PBE. While the finite el- 
ement method has the advantage of exactly representing the molecular surface (when 
appropriate mesh generation algorithms are used; see Section 5), advances in finite differ- 
ence and finite volume methods include interface discretization methods which substan- 
tially improve solution accuracy at the dielectric discontinuity surface [14, 63, 62, 20]; 
see also [54] for a similar approach using mortar elements. Boundary element meth- 
ods for the (primarily linearized) PBE are also competitive, due to algorithm advances 
for molecular surface generation (see [59, 58] and Section 5), due to emergence of fast 
multipole codes for surface integrals [34, 32, 33], and due to new techniques for nonlin- 
earity [8]. 

The remainder of the paper is organized as follows. In Section 2, we give a brief 
derivation of the standard form of the PBE, and then examine the two-scale regulariza- 
tion in [11]. We then describe a second distinct regularization and illustrate why it is 
superior to the original approach as a framework for developing numerical methods. We 
then quickly assemble the cast of basic mathematical results needed for the second reg- 
ularization, which do not immediately follow from the results established in [1 1] for the 
original regularization. In Section 3, we describe an adaptive finite element method based 
on residual-type a posteriori estimates, and summarize some basic results we need later 
for the development of a corresponding convergence theory. In Section 4, we develop 
the first AFEM contraction-type result for a class of semilinear problems that includes 
the PBE, substantially improving the AFEM convergence result given in [11]. We also 
include a discussion of our mesh generation toolchain in Section 5, which plays a key 
role in the success of the overall adaptive numerical method. Numerical experiments are 
conducted in Section 6, where stability of the regularization scheme and convergence of 
the adaptive algorithm are both explicitly demonstrated numerically, in agreement with 
the theoretical results established in the paper. We summarize our results in Section 7. 

2. The Poisson-Boltzmann Equation (PBE) 

The PBE can be derived in various ways based on the statistical description of a system 
of charged particles in electrolytes [36, 48]. A well-known derivation starts with the 
Poisson equation for the electrostatic potential <p = 4>{x) induced by a charge distribution 

P = p(x): 

-V • (eV0) = -p, 

where e = e(x) is a spatially varying dielectric constant and e is the dielectric permit- 
tivity constant of a vacuum. The charge distribution p may consist of fixed charges pf 
and mobile charges p m . The fixed charge distribution pf represents the partially charged 
atoms of the molecules immersed in the aqueous solution; the mobile charges p m models 
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the charged ions in the solution. With this perspective, the fixed charge distribution pj 
is independent of the potential <fr. The charge distribution p m of mobile ions, however, 
depends on the potential <p following the Gouy-Chapman or Debye-Hiickel theories, and 
can be modeled by a Boltzmann distribution. The two charge distributions then take the 
form 

M N 

Pm = ^2c j q j e- g ^ /kT , p f = ^2 qi S(xi) x l eQ rn . (2.1) 

3=1 i=l 

Here for p m , M is number of ion species, Cj and qj are the bulk concentration and charge 
of the j th ion, k is the Boltzmann constant and T is the absolute temperature; and for 
Pf, there are N charges located at Xi in the molecule region VL m and carrying charge % 
where 8(xi) is the delta function centered at Xi. This gives rise to the full or nonlinear 
PBE: 



. / N M \ 

- v ■ (ev^) = — y asfa) + y . a.i) 

6 ° \U 3=1 J 

A number of variations of the PBE can be derived under appropriate assumptions. For 
example, for a symmetric 1 : 1 ionic solution (two ions species with same but opposite 
charge) with M — 2, bulk concentration Cj = c and charge qj = (— l) J g for j = 1, 2, 
equation (2.2) reduces to: 

_ V • (eV0) + -2cq sinh ( = - f>5(^)- 
e o \kTJ e ^ 

We now introduce a dimensionless electrostatic potential u = qcj)/{kT), and the so- 
called Debye length l D = ^/jj^p . and define the modified Debye-Huckel parameter 

to be k = I /Id- After scaling the singular charges we can write the final form of the 
Poisson-Boltzmann equation as: 



- V • (eVu) + k 2 sinhtt = /, (2.3) 



N 



where f =} J Zi5(xi), with z { = An qqi/(e kT). 



i=l 



As analyzed in [11], since the singular function / does not belong to i/" 1 (O), equa- 
tion (2.3) does not have a solution in H 1 , or at least the equation does not have a weak 
formulation involving the H 1 as the test space. Consequently, standard numerical meth- 
ods for elliptic equations are not guaranteed to produce numerical solutions which con- 
verge to the exact solution to the PBE in the limit of mesh refinement, and numerical 
evidence suggests that in fact standard methods fail to converge. We now discuss two 
regularization schemes for the PBE which have not only been the basis for the new solu- 
tion and approximation theory results for the PBE appearing in [11], but also provide a 
robust framework for constructing provably convergent numerical algorithms. 

2.1. A Natural Regularized Formulation. The first scheme is motivated by the phys- 
ical interpretation of the solution to PBE and decomposes the solution into two compo- 
nents, based on the distinct solvent region £l s and molecular region Vt m in the model. 
This spatial decomposition of the domain Vt, as well as the interface Y between VL S and 
VL m , is depicted in Figure 1. The component of the solution, which will have singu- 
larities but will be representable in closed-form, is called the self-energy corresponding 
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Figure 1 . Illustration of the solvent region Q s , the molecular region Q m , 
the interface T, and the two distinct dielectric constants e s and e m in the 
two regions. 

to the electrostatic potential. The second component, which will be much more well- 
behaved but will not have a closed-form representation, corresponds to the screening of 
the potential due to high dielectric and mobile ions in the solution region. This natural 
decomposition provides a regularization scheme was proposed and explored numerically 
in [21, 64]. In this scheme, the singular component u s of the electrostatic potential is 
identified as the solution of the following Poisson equation, the solution of which can be 
readily assembled from the Green's functions (cf. [46]): 

N N 1 

- V ■ (e m W) = u s := V *~ (2.4) 

i=l i=l 

Subtracting (2.4) from (2.3) gives the equation for the regular component u: 

- V • (eVu) + k 2 sinh(w + u s ) = V • ((e - e m )V« s ). (2.5) 

Since k vanishes in Vt m and e — e m is nonzero only in region Vt s , the right hand side 
term V • ((e — e m )V« s ) belongs to H~ l , and a standard i^-weak formulation of (2.5) is 
well-defined. 

A variational argument can be used to show existence and uniqueness of a weak so- 
lution to (2.4) in H 1 (see [11] for this argument, and also [23] for a similar argument 
in the case of an alternative regularization). A priori L°° estimates for the solution are 
established in [11], which are critical to the development of a priori error estimates for 
Galerkin (e.g. finite element, wavelet, and spectral) approximations of the regular com- 
ponent, and are also critical to the convergence results for both uniform and adaptive 
finite element methods developed in [11]. This two-scale decomposition framework is 
at the heart of the solution theory, approximation theory, and convergence results for 
adaptive finite element methods for the PBE developed in [11]. 

2.2. An Alternative Regularized Formulation. Since the singular component repre- 
sents the Coulomb potential in the low dielectric environment, it is always much larger 
than the real potential in where the dielectric constant is high and strong ion screen- 
ing exists. As a result, the regular component is also much larger in magnitude than 
the full potential in Q s , and the decomposition in Section 2.1 can produce an unstable 
numerical scheme. More precisely, relatively small error in the numerical solution of 
regular component could lead to large relative error in the full potential, as illustrated in 
the following example. 
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Example 2.1. Let VL m be a unit ball with a unit positive charge at the origin. The di- 
electric constants are e m = 2 and e s = 80 inside and outside the ball, respectively. Let 
the modified ionic strength k = 0. This so-called Born Ion problem admits an analyti- 
cal solution for the full potential u(r) = — — h (— — ) in fl m and u(r) = — in Q s 

where r = J x 2 + y 2 + z 2 . Since the singular component u s (r) = —, it follows that 
the regular component 

u(r) = u(r) - u s (r) = —{!-—) = -39— = -39u(r). 

We assume that the singular component is computed analytically. Suppose that the nu- 
merical solution of u(r) carries a relative error e = 3%\u(r)\, and assume that this is the 
only source of numerical error. This implies the relative error of the final full potential is 
impacted as 

W 3% ^ |e| _ 0.03Kr)| _ 0.03 x 39|i2(r)| = 



\u(r)\ \u(r)\ \u(r) 

This suggests relative error of 3% in the numerical solution ofu(r) will be amplified by 
39 times in the relative error of the full potential u{r) when u{r) is added to the analytical 
solution u s {r). 

Example 2.1 indicates that unless the regular component is solved to high accuracy, 
the full potential could be of low quality, and numerical algorithms based on the decom- 
position may fail. A second decomposition scheme we now examine demonstrates more 
satisfactory numerical stability. Proposed in [14], this decomposition splits the potential 
into three parts in the molecular region only. The first component is the singular com- 
ponent u s defined by (2.4). The second component u h is the harmonic extension of the 
trace of the singular component u s on the molecular surface into the interior of the mol- 
ecule; it is completely determined by the singular component u s and the geometry of the 
molecular surface through the harmonic equation 

- Au h = in Q m , u h = -u s on T, (2.6) 

where T is the interface between Vt m and Vt s . Then we set u s +u h = in O s . By definition 
of this extension is continuous across the interface. So the complete decomposition 
reads 



, e - £ s / Sj h-, \u s + u h +u mVL m 
u = u-\ (u + u ) = < . (2.7) 



u mVt 



s 



sinh 


U 


= o, 


in Vt 


[u 


I" 


= o, 


onT, 


du 


r 




onT, 




u 


= o, 


\x\ — > 



In this decomposition, the regular component u is defined as an interface problem 
-V-(eVu) 

ktL = 0, on T, 

witng r .— e m — ^ — | r , 

where n r is the unit out normal of the interface T, and [•] denotes the jump of enclosed 
quantity on the given interface as [t>] r = lim^o v(x + tn r ) — v(x — tn r ). The second 
interface condition (2.8) arises from continuity of flux in (2.3). The singular component 
u s is given by (2.4), whereas computing u h is trivial using finite element or boundary 
integral methods. Therefore, we assume u s and u h are known in the following discussion, 
and are smooth on Y. Since the singular component u s is only applied in the interior of 
the molecule region, a discontinuity appears in the remaining component of the potential 
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on the molecular surface. The harmonic component u is introduced to compensate for 
this discontinuity using harmonic extension, so that the regular component as defined by 
equation (2.8) is continuous on the molecular surface. 

Since no decomposition of the potential occurs in Vt s , error in numerical solutions of 
u are not amplified in the full potential. While mathematically equivalent to the decom- 
position in [11], this alternative three term-based splitting regularization is potentially 
numerically more favorable than the original decomposition. The implementations of 
this scheme using finite difference interface methods [14, 20] have proven that it can 
significantly improve the accuracy of the full potential. Mirroring the general plan taken 
in [11], we will use this alternative decomposition as the basis for an analysis of the 
regularized problem, for the development of an approximation theory, and for the devel- 
opment of a practical, provably convergent adaptive method. 

A final difficultly in solving the regularized form of the PBE in (2.8) (and other forms 
of the PBE) is that the computational domain is all of space. It is standard to truncate 
space to a bounded Lipschitz domain Vt by posing some artificial (but highly accurate) 
boundary condition on dVt. For simplicity, one chooses Vt to be a ball or cube contain- 
ing the molecule region. The solvent region is then defined as Vt s fl f2, which will also 
be denoted by Vt s without the danger of confusion. There are various approaches to the 
choosing boundary condition on dVt; using the condition u — g is standard, where g can 
be obtained from a known analytical solution to some simplification of the linearized 
PBE, and can be chosen to be a smooth function on the boundary. Far from the molecule 
region, such analytical solutions provide a highly accurate boundary condition approxi- 
mation for the PBE on the truncated domain. For other possible constructions of g, see 
[11,7, 23] and the references cited therein. Finally, we end up with the regularized PBE 
(or RPBE) in a bounded domain fi, which becomes the focus for the remainder of the 
paper: 



Our main goals for the remainder of the paper are to: 

1) Establish a priori L°° estimates for (2.9), leading to a standard argument for well- 
posedness of the continuous and discrete problems. Most other mathematical re- 
sults in the article hinge critically on these a priori estimates. 

2) Develop a general approximation theory for (2.9) by establishing a priori error 
estimates for Galerkin methods, giving convergence of finite element and other 
methods; 

3) Develop a practical adaptive finite element method for (2.9) and prove that it is 
convergent; 

4) Develop practical mesh generation algorithms for the domains arising in (2.9) that 
meet the needs of our finite element methods. 

We note that there are two distinct interface conditions in (2.9), which appears to give it 
an unusual formulation. However, the first interface condition [u] = will be automati- 
cally satisfied by standard constructions of C° finite element spaces. The second interface 
condition will be embedded into the weak form of equation (2.9) in a natural way, so that 
in fact both interface conditions are quite easily and naturally incorporated into finite 
element (as well as wavelet and spectral) discretizations. Although fairly complicated 
schemes arise when considering the regularization approach with finite difference and 




(2.9) 
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finite volume methods, the interface conditions can be enforced with these discretization 
as well (cf. [52, 39, 30, 29]). 

2.3. A priori L°° -Estimates and Well-posedness. A priori L°° estimates for the solu- 
tion u to the regularized PBE (2.9) are the critical component of the key mathematical re- 
sults we need to have in place for the development of a reliable adaptive method, namely: 
(1) well-posedness of the continuous and discrete regularized problems; (2) a priori er- 
ror estimates for Galerkin approximations; (3) a posteriori error estimates for Galerkin 
approximations; and (4) auxiliary results for establishing convergence (contraction) of 
AFEM. The regularized equation (2.9) governing u derived in Section 2.2 differs signif- 
icantly from the decomposition used in [11], and as a consequence we now derive the a 
priori L°° estimates. 

In what follows, we use standard notation for the L P (G) spaces, 1 ^ p ^ oo, with 
the norm || • \\ Pi q on any subset G C W 1 ; we use standard notation for Sobolev norms 
IMkp.G — ll M llvF fc >p(G) where the natural setting here will be p = 2 and k = or k = 1. 
For any functions v e L P {G) and w £ L q {G) for p, q ^ 1 with - + - — 1, we denote the 

pairing (v, w)c as (v, w)q '■— J G vwdx. If G = Vl then we also omit it from the norms 
(or pairing) to simplify the presentation. 

To begin, define an affine subset of H 1 ^) as H^Q) := { v e H 1 ^) : v = g on d£l }, 
and then define X g := { v G H^(Q) : e v ,e' v e L°°(fi) }, with X denoting the case 
when g = 0. A weak formulation of equation (2.9) reads: Find u G X g such that 

a(u,v) + (b(u),v) = (g T ,v) T , Vt; e H^Q), (2.10) 

where a(u,v) = (eVu,Vt>), (b(u),v) = (k 2 sinh(u), v ) and (gr,v)r = J r grvds. It is 
easy to verify that the bilinear form in (2.10) satisfies: 

m\\u\\l 2 ^ a ( u , u ) i a {u, v) ^ M||"u||i ) 2||i>||i,2, Vw, v e H (fl) , (2.11) 

where < m ^ M < oo are constants depending only on the maximal and minimal 
values of the dielectric and on the domain. The properties (2.11) imply the norm on 
Hq(Q) is equivalent to the energy norm ||| • ||| : Hq(CI) — > E, 

|||m||| 2 = a(u,u), m||n||f 2 < |||m||| 2 ^ M||m||^ 2 . (2.12) 

To establish a priori L°° estimates, we further split the solution u to (2.9) into solu- 
tions of two sub-problems. The first sub-problem is a linear elliptic interface problem; 
estimates on solutions to this problem are then utilized in the analyzing the second sub- 
problem, which is a nonlinear elliptic problem without interface conditions. The second 
sub-problem is then analyzed using a cut-off function argument that exploits a weak for- 
mulation of the maximum principle. More precisely, let u = u l + u n , where u l E Hg(St) 
satisfies the linear elliptic equation 

a(u l ,v) = (g r ,v) r , Wv e H^Q), (2.13) 

and u n E Xq satisfies the nonlinear elliptic equation: 

a{u n ,v) + (b(u l + u n ),v) = 0, VveH^n), (2.14) 

where we note that the sum u = u l + u n is then the desired solution to the RPBE (2.10). 
It is easy to see that the linear part u l is the solution to the interface problem: 

-V ■ (eW) = in Q 

u l \ dn = g, and 



, 8u l 
' 9rtp 
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while the nonlinear part u 11 is the solution to the (homogeneous) semilinear equation 

f -V ■ (eVw n ) + k 2 sinh(w n + u l ) = in Q, 
\ u n = ondfl. 

Existence and uniqueness of u l solving (2.13) follows by standard arguments; further- 
more, if the interface T to be sufficiently smooth (e.g. T is C 2 ), then u l G follows 
immediately from known regularity results for linear interface problems (cf. [11, 4, 9, 
12, 43]). This makes possible a priori L°° estimates for the nonlinear component, and 
subsequently the entire regularized solution. To this end, define 

a = argmax ( k 2 sinh(c + sup u{) ^ ) , a n = min(a', 0), (2.15) 
ft = argmin (k 2 sinh(c + inf U j) ^ oY (3 n = max(/3', 0). (2.16) 

Lemma 2.2 (A Priori L°° Estimates). Suppose that the solution u l to (2.13) satisfies 
u l G L°°(Q), and let u n be any weak solution of (2.14). Ifa n , (3 n G K are as defined 
in (2.15)-(2.16), then 

a n ^u n ^ (3 n , a.e. inVL. (2.17) 

Proof. The short proof is similar to that in [11, 46], which we include for completeness, 
due to its critical role in the results throughout the article. We first define 

= ( u n - P n ) + = max(u n - P n , 0), = {u n - a n )~ = max(a" - u n , 0). 

Since f3 n ^ and a n ^ 0, it follows (cf. [46]) that 0,0 G H^(Q), and can be used as 
pointwise non-negative (almost everywhere) test functions. For either = or = —0, 
we have 

(eVu n , V0) + (k 2 sinh(u n + u l ), 0) = 0. 
Note ^ in Q and its support set is y = {x G | -u n (x) ^ /3 n }. On 3^, we have 

sinh(u n + u') ^ k 2 sinh(/3' + inf u l ) ^ 0. 

Similarly, —0 ^ in with support ^ = {i 6 H | u n (x) ^ a n }. On 3^, we have 
k 2 s'mh(u n + u l ) ^ k 2 sinh(a' + sup u l ) ^ 0. 

x£fl s 

Together this implies both 

^ (eW\V0) = (eV(« n -/3 n ),V0) ^ (inf e(x))\\ V4>\\\ ^ 0, 

ISO 

^ (eW\V(-0)) 2 = (eV(« n -M n ),V0) 2 ^ (inf £(x))||V0|| 2 ^ 0. 

— — x£Q — 

Using the Poincare inequality we have finally ^ ||0||i,2 < || V0||2 ^ 0, giving = 0, 
for either = or = —0. Thus a n ^ u n ^ (3 n in fi. □ 

We have therefore shown that any solution u G H 1 (fl) to the regularized problem (2.10) 
must lie in the set 

[a, /3] lj2 := { u G H\n) : a ^ u < /3 } C ff 1 ^), 
where a, /? G M. are 

a = a n + inf ui, (3 = (3 n + sup u\. 
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Since this ensures u G L°°(Cl), which subsequently ensures e u G L 2 (Q), we can replace 
the set X g with the following function space as the set to search for solutions to the 
RPBE: 

V = { u G [a, /3]i >2 : u = g on dtt } C X g C H\Q). 
Our weak formulation of the RPBE now reads: 

Findw G V such that a(u,v) + (b(u),v) = (gr,v) r , Vf G Hq(Q). (2.18) 

Note that in general V is not a subspace of H 1 ^) since it is not a linear space, due 
to the inhomogeneous boundary condition requirement. However, as remarked above, 
standard results for linear interface problems imply existence, uniqueness, and a priori 
L°° bounds for u l solving (2.13), leaving only the equation (2.14) for the remainder u 11 . 
Therefore, (2.18) is mathematically equivalent to 

Find u n G U such that a(u n ,v) + (b(u n + u l ),v) = Vt> G Hq(Q), (2.19) 

where 

U = { u G Hl(£l) : u_ < u < u + } C F^O), 

with m_ = a n and w + = /3 n from Lemma 2.2. We now have a formulation (2.19) that 
involves looking for a solution in a well-defined subspace U of the (ordered) Banach 
space X = Hq(Q), and are now prepared to establish existence (and uniqueness) of the 
solution. The argument we use below differs significantly from that used in [11, 23] for 
the original regularization. 

Theorem 2.3 (Existence and Uniqueness of Solutions to RPBE). Let the solution u l to 
(2.13) satisfy u l G L°°(Q). Then there exists a unique weak solution u n G U to (2.19), 
and subsequently there exists a unique weak solution u G V to the RPBE (2.18). 

Proof. We follow the approach in [1 1, 23]. We begin by defining J : U C H] (Q) -)> E: 

J(u) = I -\ Vm| 2 + k 2 cosh(-u + u l ) dx. 
Jn 2 

It is straight-forward to show that if u is the solution of the optimization problem 

J{u) = inf J(v) < J(v), Vv G U, (2.20) 

V&J 

then u is the solution of (2.19). We assemble some quickfacts about Hq(VL), U C Hq(Q), 
and J. 

1) Hq(Q) is a reflexive Banach space. 

2) [/ is nonempty, convex, and topologically closed as a subset of Hq(Q,). 

3) J is convex on U: J(Xu+(l-X)v) < AJ(m) + (1-A)J(u),Vu,t; G 17, A G (0,1). 

By standard results in the calculus of variations (cf. [46]), we have existence of a solution 
to (2.20), and hence to (2.19) and (2.18), if we can establish two additional properties of 
J: 

4) J is lower semi-continuous on U : J(u) ^ liminfj^oo J(uj), Vnj — >■ u G Z7. 

5) J is coercive on Z7 : J(«) ^ Co||u||f 2 — Ci» e U ■ 

That J is lower semi-continuous (and in fact, has the stronger property of weak lower 
semi-continuity), holds since J is both convex and Gateaux-differentiable on U (cf. [46] 
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for this and similar results). That J is coercive follows from cosh x ^ and the Poincare 
inequality 

f S Si Ob i 

J(u) = / -|Vw| 2 + k 2 cosh(w + u l )dx ^ inf -^-\u\\ 2 

^ inf £ M ( I| u |2 + * ^ C'ollitHfa, 

^ xen 2 V 2 V J 

with C = (inf xg n e(x)) ■ min{l/4, l/(4p 2 )}, where p > is the Poincare constant. 
It remains to show u is unique. Assume there are two solutions u\ and u 2 . Subtract- 
ing (2.19) for each gives a{u\ — u 2 ,v) + (b(ui + u l ) — b(u 2 + u l ),v) = Vt> G 
//o(f2). Now take t> = U\ — u 2 ; monotonicity of the nonlinearity defining b ensures that 
{b(u\ + u l ) - b(u 2 + u l ),ui - u 2 ) ^ 0, giving ^ a{u\ - u 2 ,ui - u 2 ) ^ 2C ||wi - 
w 2 1| 1,2 ^ 0, where Co is as above. This can only hold if u\ = u 2 . □ 

In summary, there exists a unique solution u E V C to the RPBE prob- 

lem (2.18), with compatible barriers ?i_ and u + E L°° satisfying 

— oo < U- ^ u ^ u + < oo, a.e. in Q. 

Moreover, these pointwise bounds combined with a Taylor expansion give that for any 
u,w E V and any v E Hq(Q), the nonlinearity satisfies a Lipschitz condition: 

(b{u) -b{w),v) ^ K\\u-w\\ 2 \\v\\ 2 , (2.21) 

where K = sup xg r U )U i \\k 2 cosh(x)||oo < oo is a constant depending only on the do- 
main, the ionic strength of the solvent (embedded in the constant k), and other physical 
parameters. 

3. Finite Element Methods (FEM) 

In this section, we consider (mainly adaptive) finite element methods for the regular- 
ized problem (2.9). For simplicity, we assume VL be a bounded polygon domain, and we 
triangulate Vt with a shape regular conforming mesh Th- Here h = /i max represents the 
mesh size which is the maximum diameter of elements in Th- We further assume that 

Assumption: Al the discrete interface Th approximates the original interface Y to the 
second order, i.e., d(T,Th) ^ ch 2 . 

The mesh generator discussed in Section 5 provides a practical tool for generating meshes 
with this type of approximation quality for the interface. Given such a triangulation Th, 
we construct the linear finite element space 

V(%) := {v E H l (0) : v\ T E V X {t), Vr E %}■ 

Since we may choose g to be a smooth function on dfl, the Trace Theorem (cf. [1, 
46]) ensures there exists a fixed function ud € H l {VL) such that «o = g on dfl in 
the trace sense. Let H^Q) := Hq(Q) + w fl be the affine space with the specified 
boundary condition, and Vd(%,) '■= V(%) H H^Q) be the finite element affine space of 
In particular, we denote V (Th) '■= V(Th) H Hq{Q). For simplicity, we assume 
the boundary condition g can be represented by exactly. In practical implementation, 
we will construct an interpolant of ud having sufficient approximation quality such that 
using the interpolant in place of ud will not impact the order of accuracy of the algorithm 
we build below for approximating the solution u to the regularized problem (2.9). A 
Galerkin finite element approximation of (2.10) takes the form: Find uu E VniTh) such 
that 

a(u h ,v) + (b(uh),v) = (gr,v), Wv E V (%). (3.1) 
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The primary concerns for this type of approximation technique are the following four 
mathematical questions regarding the Galerkin approximation u h : 

1) Does Uh satisfy discrete a priori bounds in L°° and other norms, so that the non- 
linearity can be controlled for error analysis? 

2) Does Uh satisfy quasi-optimal priori error estimates, so that the finite element 
method will converge under uniform mesh refinement? 

3) Can one ensure AFEM (non-uniform mesh refinement) convergence lim^oo Uk = 
u, where u^ is the Galerkin approximation of u at step k of AFEM? 

4) Can one produce u h at each step of the (uniform or adaptive) refinement algorithm 
using algorithms which have optimal (linear) or nearly optimal space and time 
complexity? 

The first two questions were answered affirmatively in [11] for the first regularized for- 
mulation. We give only a brief outline below as to how the arguments for answering 
the first two questions can be modified to establish the analogous results for the sec- 
ond regularized formulation here. The third question was partially answered in [11] for 
the first regularization, but we give an improved, more complete answer to this ques- 
tion below and in Section 4, which is one of the main contributions of the paper. Re- 
garding the fourth question, due to the discontinuities in the dielectric and the modified 
Debye-Huckel parameter, one must take care to solve the resulting nonlinear algebraic 
systems using robust inexact global Newton methods, combined with modern algebraic 
multilevel-based fast linear solvers in order to produce an overall numerical solution 
algorithm which is reliable and has low-complexity. We do not consider this question 
further here; see [26, 3, 2, 23] for a complete discussion in the specific case of the PBE. 

3 .1. Discrete L°° Estimates and Quasi- Optimal A Priori Error Estimates. For com- 
pleteness, we quickly answer the first two questions by stating a result, giving only a 
very brief outline of how the result is established for the new regularization, based on 
modifying the analogous arguments in [11]. We then focus entirely on the new AFEM 
contraction results which require a more complete discussion. To state the theorem, the 
following assumption is needed. 

Assumption: A2 For any two adjacent nodes i and j, assume that 

ciij = (eV0i, V<f>j) < l r l> with c > °> 

eijCr 

where e^- is the edge associated with these nodes, (pi and <pj are the basis functions 
corresponding to nodes i,j respectively, and |r| is the volume ofreT. 

Theorem 3.1. If Assumption A2 holds and h is sufficiently small, the solution to (3.1) 
satisfies: 

Halloo ^ c, 

where C is independent ofh. Moreover, the quasi-optimal a priori error estimate holds: 

\\u-u h \\i >2 < inf ||m-u||i,2- 

veV D {T h ) 

Proof. The proofs of both inequalities are similar to the proofs of the corresponding 
results in [11], with adjustment to handle the new stabilized splitting. The second result 
hinges critically on the Lipschitz property (2.21), which in turn relies on the first result 
together with the continuous L°° estimates established in Lemma 2.2. □ 
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3.2. Adaptive Finite Element Methods (AFEM). Adaptive Finite Element Methods 
(AFEM) build approximation spaces adaptively; this is done in an effort to use nonlinear 
approximation so as to meet a target quality using spaces having (close to) minimal 
dimension. AFEM algorithms are based on an iteration of the form: 

SOLVE -> ESTIMATE -> MARK -> REFINE 

which attempts to equi-distribute error over simplices using subdivision driven by a pos- 
teriori error estimates. Given an initial triangulation 7o, and a parameter 9 E (0, 1], our 
particular AFEM generates a sequence of nested conforming triangulations Th, h > 0, 
driven by some local error indicator r](uh, t), which gives rise to a global error indicator 
v( u hi Th) - Schematically, the adaptive algorithm consists of a loop of the following main 
steps: 

1) u h := SOLVER) . 

2) {v(uh,r)} T er h ■= ESTIMATE^, %) . 

3) M h := MARK({n(u h ,T)}rer h ,Th,e). 

4) % :=REF\NE(T h ,M h ,£). 

In practice, a stopping criteria is placed in Step (2) to terminate the loop. 
We will handle each of the four steps as follows: 

1) SOLVE: We use standard inexact Newton + multilevel to produce U E Vd(T) 
on triangulation T (cf. [24]). To simplify the analysis here, we assume that the 
discrete solution U is calculated exactly (no round-off error). Given a triangulation 
T, this defines the procedure: 

U := SOLVE(T). 

2) ESTIMATE: Given a triangulation T and a function U E Vd(T), we compute the 
elementwise residual error indicator: 

MU,r)} Ter := ESTIMATE^, T). 

3) MARK: We use the standard "Dorfler marking": 

Given 9 E (0, 1], we construct a marked subset of elements 

M := MARK({r)(U,T)} TeT ,T,6) C T, 

such that: 

ij(U,M)>9ri(U,T). (3.2) 

The residual-type error indicator r)(U,At) over sub-partition M C T will be de- 
fined precisely in Section 3.3, 

4) REFINE: We use standard non-degenerate bisection-to-conformity methods with 
known complexity bounds on conformity preservation (cf. [47]). In particular, 
given a triangulation T, and marked subset M. C T, and an integer i ^ 1, we 
produce 

X := REFINE(T, AV), 
a conforming refinement of T with each simplex in M. refined at least t times. 

3.3. Residual A Posteriori Estimates. To make precise the residual-type indicator, we 
first introduce some standard notation for the relevant mathematical quantities, and then 
employ and establish some of its properties. 
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7o = Initial conforming simplex triangulation of C W 1 . 

Th = Conforming refinement of Tzj at the previous step of AFEM. 

h T = The diameter of simplex r G T. 

= The normal vector to face FofreT. 

u T = U { f G T | rf|T 7^ 0, where r G T}. 

Wf = |J { t G T | Ff]f ^ 0, where F is a face of r G T }. 
We can now define the following error indicators: 

V 2 {u h ,r) := hl\\b{u h )\\l T + ^^ h F \\n F -[eVu h }\\l F + ^ MMIjyr, 

Fc9r Fcdmr h 

(3.3) 

osc 2 (« ft ,r) := /i*||V« A ||2 jT + ^2 h F\\9r - gr\\l,Fi ( 3 - 4 ) 

where g r is the piecewise average on each face F C T h . For any subset S C T, the 
cumulative indicators are defined as: 

r? 2 ^, S) := r/ 2 ^, r), osc 2 ^, 5) := osc 2 (w fe , r). 

From these definitions follows the monotonicity properties: 

v (v,%) < ^,T), We^T), (3.5) 

osc(w,7;) < osc(v,T), Vf G Vd(T), (3.6) 

for any refinement 7^ of T. We have then the following global upper bound from [24, 50]; 
the lower-bound is also standard and can be found in e.g. [50]. 

Lemma 3.2 (Upper and lower bounds). Let u and Uh be the solutions to (2. 10) and (3.1), 
respectively. If the mesh conditions in Assumptions Al and A2 hold, then there exists 
constants C\ and C 2 , depending only on % and the ellipticity constant, such that the 
following global upper and lower bounds hold: 

\\\u-u h \\\ 2 < C x r] 2 {u h ,Th), (3.7) 
C 2 r] 2 (u h ,r) < \\\u - u h \\\l T + osc 2 (u h ,u T ), (3.8) 
where III fill 2 , = f e\Vv\ 2 dx. 

I" IHWt Juj t I I 

Proof. Similar to [11, Theorem 7.1], the proof follows the idea of [50] by noticing 



\\b(u h ) - b{u h )h,r ^ \\K 2 (smh(u h ) - sinh(w^))|| 2 ,T 

2 



< Slip \\K, COsh(x)\\oo,r\\ U h - Uh\\2,T 

XG[n_,w + ] 

< C(K,U-,U + )h T \\Vu h \\ 2 ,r, 

where we have used Theorem 3.1. The remaining proof is the same as [11, Theorem 
7.1]. " □ 

Note for the convergence analysis here, we do not need the lower bound (3.8). 

4. Convergence of AFEM 

We now develop a convergence analysis of the AFEM iteration by showing contrac- 
tion. We must establish two additional key auxiliary results first: an indicator reduction 
result, and a quasi-orthogonality result, which generalize two analogous results for the 
linear case in [10] to a class of nonlinear problems that includes the Poisson-Boltzmann 
equation. 
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4.1. An Indicator Reduction Lemma. Here we establish a nonlinear generalization of 
the indicator reduction result from [ 1 0, Corollary 4.4] . First we prove a local perturbation 
result for the nonlinear equation (cf. [10, Proposition 4.3]). We then establish an indicator 
reduction result. 

Let us first introduce a type of nonlinear PDE-specific indicator: 

V 2 (D,t) := ||e||^ + ^ sup \\k 2 cosh( X ) ||L, T 

XE[u-,u+] 



For any subset S C T, let ^(D, <S) := max Te s{r](D, r)}. By the definition, it is obvious 
that ??(D, T) is monotone decreasing, i.e., 

V(D,%) ^V(D,T) (4.1) 

for any refinement T* of T. 

We now establish (see also [27]) a nonlinear generalization of the local perturbation 
result appearing as in [10, Proposition 4.3]. This is the key result in generalizing the 
contraction results in [10, Proposition 4.3] to the semilinear case. 

Lemma 4.1 (Nonlinear Local Perturbation). Let T be a conforming partition satisfies 
Assumptions AI and A2. For all r e T and for any pair of discrete functions v , w G 
[u-, u+) PI Vd(T), it holds that 

i](v,t) ^ r}(w, t) +A 1 »7(D,t)||t; - 1,2,1*. , ( 4 - 2 ) 

where Ai > depends only on the shape-regularity of To, and the maximal values that b 
can obtain on the L°° -interval [u-, u + ]. 

Proof. By the definition (3.3) of 77, we have 

t)(v,t) < r)(w,T) + h T \\b(v)-b(w)\\ 2}T + - hMn F -[eV(v-w)]\\ 2 , F 



2 

Fcdr 



Notice that 



||fo(w)-6(w)|| 2 , T = ||K 2 (smh(v)-sinri(u;))||2, T ^ sup \\k 2 cosh(x)||oo,r ) ||f-w||2,r- 
On the other hand, we also have 

_ 1 

\\n F ■ [eV(f - w))\\ 2>F ^ Helloo^/ir 2 ||Vv - Vw\\ 2>u}t . 
Therefore ,we get the desired estimate for rj. □ 

Based on Lemma 4.1, we have the following main estimator reduction (see also [27]), 
which generalizes the linear case appearing in [10, Corollary 4.4]. 

Lemma 4.2 (Nonlinear Estimator Reduction). Let T be a partition which satisfies the 
mesh conditions in Assumptions A 1 andA2, and let the parameters 9 G (0, 1] and I ^ 1 
be given. Let M = MARK({r](v,r)} TeT ,r,9), and let % = REFINE(T, M, £). If 
Ai — (d+ l)A.j/£ with Ax from Lemma 4.1 and A = 1 - 2~ (£/d) > 0, then for all 
v G [u-, u + ] n Vd(T), f* G [u~, u+] H Vjj(%), and any 5 > 0, it holds that 

rj 2 (v*, %) < (1 + 5)[ V 2 (v, T) - \ V 2 (v, M)} + (1 + r^A^^D, 73) Ik - vf. 

Proof. We follow the proof in [10, Corollary 4.4] closely. We first apply Lemma 4. 1 with 
v and t>* taken to be in Vd{%). This gives 

V(v*,n) 2 «C {l + 5)r}{v,n) + (1 + r 1 )A^(D,r*)||^ - v\\ 1)2l ^ Vr» G %, 
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after applying Young's inequality with 5 > 0. We now sum over the elements r* G %, 
using the fact that for shape regular partitions there is a small finite number of elements 
in the overlaps of the patches co Tt . This gives 

V(v*, %f < (1 + S)r](v, %) + (1 + r 1 )A^(D, %) \\\v* - v\\, 

where we have also used the equivalence (2.12). 

Now let v G [w_,it+] fl Vd(7"), a short argument from the proof of Corollary 4.4 
in [10] gives 

r] 2 (v, %) < if(v, T\M) + 2-^ V 2 (v, M) = rf{v, T) - \rf{v, M). (4.3) 

Finally, the monotonicity properties r)(D,%) ^ r](D,%), combined with (4.3) yields 
the result. □ 

4.2. Quasi- Orthogonality for Nonlinear Problems. Following [27], we now establish 
a quasi-orthogonality result that represents the last technical result needed to generalize 
the convergence framework from [10] to the nonlinear case. 

Lemma 4.3 (Quasi-orthogonality). Let u be the exact solution to equation (2.10), and Uh 
be the solution to (3.1) on a partition Th which satisfies the conditions in Assumptions Al 
and A2. Assume that there exist a <Jh > with ah — > as h — > such that 

\\u - u h \\ 2 < <Th\\Vu - Vu h \\ 2 , (4.4) 

Then there exists a constant C* > 0, such that for sufficiently small h, we have 

\\\u - u h \\\ 2 < A h \\\u - u H \\\ 2 - \\\u h - u H \\\ 2 , (4.5) 

where A h — (1 - C'd,^) -1 > with K = sup xg[u _ u+] \\k 2 cosh(x) ||oo- 

Proof. We compute the energy norm: 

— uh\\\ 2 = a{u — uh, u — Uh) = a(u — Uh + Uh — Uh, u — Uh + u h — u h) 

= \\\u - u h \\\ 2 + \\u h - u H \\\ 2 + 2a(u - u h , u h - u H ). 

By the definition of u and u h , we have 

a(u - u h , v h ) + (b(u) - b(u h ),v h ) = 0, Vv h G V (T h )- 

In particular, this holds for Vh = Uh — %■ By this relation, we obtain 

I it - u H \\\ 2 = \\\u - u h \\\ 2 + 1 u ft - u H \\\ 2 + 2(K 2 (sinh(w) - smh(u h )),u h - u H )- 

Therefore, by Young's inequality and the assumption (4.4) we have 

|2(6(u) - b(u h ),u h - u H )\ < 2 sup ||/t 2 cosh(x)||oc||w-Wh.||2|K-«ff||2 

xe [u-,u+] 

K 2 

< 5\\u - u h \\\ + -£-\\ u h ~ u H \\\ 

K 2 

< 5C(e)cr 2 h \\\u - u h \ 2 + ^IIK - u h\\ 2 , 

for 5 > to be chosen later, and m the coercivity constant. For a h sufficient small, we 
have 

(1 - 5C7(e)a 2 )||n - u,J 2 ^ - u H f - (l - ^- J \\\u h - u H f . 

Define now the constant C* = C(e)/y/m and take 5 = Kf (^/C(e) y/mah)- We assume 
ah is sufficiently small so that 5C(e)a 2 L = C*ahK < 1. This gives (4.5) with = 
(l-C*a h K)-\ ' □ 
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We note (4.4) can be established by "Nitsche trick" under regularity assumptions (cf. 
[37]). 

4.3. The Main Convergence Result for AFEM. To establish this result, we will fol- 
low [27] and use a combination of the frameworks in [10, 38] rather than from [11]. 
This is because these frameworks are the first to handle dependence of the oscillation on 
the discrete solution itself. The quasi-orthogonality result is explicit in [38], but some- 
what hidden in [11]. The framework in [10] uses only orthogonality rather than quasi- 
orthogonality, but has a number of improvements over [38] and [11] in several respects, 
including a one-pass algorithm using only the residual indicator. 

The previous sections focused on establishing some supporting results involving a 
nesting of three spaces Xh C Xh C X, where these were abstract spaces in some 
cases, or specific finite element subspaces of H 1 . In what follows, we now consider the 
asymptotic sequence of finite element spaces produced by the AFEM algorithm, and will 
use the results of the previous sections with the subscript h in Uh and other quantities 
replaced by an integer k representing the current subspace generated at step k of AFEM. 
To simplify the presentation further, we also denote 

e k = \\\u - u k \\, E k = \\\u k - u k+ i\\\, 

Vk = v(u k ,T k ), r]k{Mk) = i](u k ,M k ), 770(D) = r] (D,T ), 

where D represents the set of problem coefficients and nonlinearity. We also denote 
V k := Vd(%) for simplicity. The supporting results we need have been established 
in §4.1 and §4.2. 

Theorem 4.4 (Contraction). Let {T k , V k , u k } k ^o be the sequence of finite element meshes, 
spaces, and solutions, respectively, produced by AFEM(9,£) with marking parameter 
9 G (0, 1] and bisection level £ ^ 1. Let T k satisfy the conditions in Assumptions Al 
and A2, and h,Q be sufficiently fine so that Lemma 4.3 holds for {T k , V k ,u k } k ^o- Then, 
there exist constants 7 > and a G (0, 1), depending only on 9, £, and the shape- 
regularity of the initial triangulation To, such that 

\\u - Mfc+if + -fvl+i «s " 2 (Ih - u kf + ivl) ■ 

Proof. We combine the frameworks in [10, 38] using the quasi-orthogonality result in 
Lemma 4.3 rather than the approach in [11] for nonlinearities. Our notation follows 
closely [10]. The proof requires the following tools: 

1) Dorfler marking property given in equation (3.2). 

2) The global upper-bound in Lemma 3.2. 

3) Estimator reduction Lemma 4.2. 

4) Quasi-orthogonality Lemma 4.3. 

In addition, some results above used indicator monotonicity properties (3.5)-(3.6) and 
monotonicity of data T] k (D). Starting with the quasi-orthogonality result in Lemma 4.3 
we have 

4+1 ^ A fe+i e fc - E l, 

which gives 

e 2 k+1 + -yr)l +1 $C A k+ ie 2 k - E 2 k + 7 ^ +1 . 
Employing now Lemma 4.2 for some 5 > to be specified later we have 

4+i + 7^+i < Afc+ie 2 : - E 2 k + (1 + 5H4 - X V l(M k )} + (1 + <T ^Ai^D)/^ 2 , 
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where A G (0, 1) as defined in Lemma 4.2. Take now 5 > sufficiently small so that we 
can ensure 7 < 1 by setting: 

o< ^w= (1+ *)Lg(D) <L (4 - 6) 

Using (4.6) in the last term leads to 

e 2 k+1 + 7^ + i < Ak+iel + (1 + 8)^1 - (1 + 5)X 1V 2 k (M k ). 
We now use the marking strategy in equation (3.2) to give 

e 2 k+1 + W k+ i < A^e 2 + (1 + 5) 1V 2 k - (1 + 5)A 7 fl 2 r/ 2 . (4.7) 

To allow for simultaneous reduction of the error and indicator, we follow [10] and split 
the last term into two parts using an arbitrary f3 G (0, 1): 



+ lv l + i < A fc+1 e 2 + (1 + 5) 1V l - 0(1 + 8)X 7 e 2 V l - (1 -/?)(! + 5)X 1 9 2 r, 2 k . 



We now the first and third terms using the upper bound from Lemma 3.7 and the expres- 
sion for 7 in (4.6), and combine the second and fourth terms as well, giving: 

4 +1 + ivl+i < (Aw - gl A^g(D) ) el + {1 + 6){1 ~ (1 ~ fl^W* 



This can be written in the form 



e 



i+i + Wk+i < <*l{6, (3)e 2 k + 7^(5, 0)r]l, 



where 



ai(5,P) = A h -5 



(3X6 



2 



CiAi% 2 (D) 



a 2 2 (5,P) = (l + 8)(l-(l-(3)Xe 



2\ 



By Lemma 4.3, we have A fc+i = (1 — C*a k+ iK) with a k+ i := <J/ lfc+1 , so that 



1-C*a fc+1 K LCiAi7?g(D). 

By the assumptions in Lemma 4.3, we can take the initial mesh so that a k+ i ^ is 
as small as we desire, or that A fc+1 is as close to one as we desire. Therefore, we can 
simultaneously pick cr fe+1 > and 5 > sufficiently small so that a\ < 1. Either this 
choice of 5 > ensures a 2 < 1 as well, or we further reduce 5 so that: 

a 2 = maxima 2 .} < 1- 
This completes the proof. □ 

5. Feature-Preserving Mesh Generation for Biomolecules 

Mesh generation from a molecule is one of the important components in finite element 
modeling of a biomolecular system. There are two primary ways of constructing molec- 
ular surfaces: one is based on the 'hard sphere' model [41] and the other is based on the 
level set of a 'soft' Gaussian function [22]. In the first model, a molecule is treated as a 
collection of 'hard' spheres with different radii, from which three types of surfaces can 
be extracted: van der Waals surface, solvent accessible surface, and solvent excluded 
surface [15, 22, 28, 41]. The molecular surface can be represented analytically by a 
list of seamless spherical patches [15, 49] and triangular meshes can be generated using 
such tools as MSMS [42]. In contrast, the 'soft' model treats each atom as a Gaussian- 
like smoothly decaying scalar function in IR 3 [6, 17, 22]. The molecular surfaces are then 
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approximated by appropriate level sets (or iso-surfaces) of the total of the Gaussian func- 
tions [6, 17]. Because of its generality, robustness, and capability of producing smooth 
surfaces, we will utilize the 'soft' model (or level set method) in our molecular mesh 
generation. 

We now briefly outline the algorithms of constructing triangular and tetrahedral meshes 
from a molecule that is given by a list of centers and radii for atoms (e.g., PQR files [16] 
or PDB files with radii defined by users [5]). More details can be found in our earlier 
work [58]. Figure 2 shows the pipeline of our mesh generation toolchain. Note that our 
tool can also take as input an arbitrary 3D scalar volume or a triangulated surface mesh 
that has very low quality. 



Input 

OR 
Input 

OR 
Input 



Atom Coordinates 



3D Volumes 



Initial Surface Meshes 



Quality Improvement 



Mesh Coarsening 



Tetra Generation 



Output 



Figure 2. Illustration of our mesh generation toolchain. The inputs can 
be a list of atoms (with centers and radii), a 3D scalar volume, or a user- 
defined surface mesh. The latter two can be thought of as subroutines of 
the first one. 



5.1. Molecular Surface Generation. In our mesh generation toolchain, a molecular 
surface mesh is defined by a level set of the Gaussian kernel function computed from a 
list of atoms (represented by centers Cj and radii r^) in a molecule as follows [6, 22, 60]: 

F(x) = 5> * \ (5.1) 

i=l 

where the negative parameter Bi is called the blobbyness that controls the spread of char- 
acteristic function of each atom. The blobbyness is treated in our work as a constant pa- 
rameter (denoted by B ) for all atoms. Our experiments on a number of molecules show 
that the blobbyness at —0.5 produces a good approximation for molecular simulations. 

Given the volumetric function -F(x), the surface (triangular) mesh is constructed us- 
ing the marching cube method [31]. Figure 3(A) shows an example of the isosurface 
extracted using this method. From this example, we can see that: (a) the isosurfac- 
ing technique can extract very smooth surfaces, but (b) many triangles are extremely 
"sharp", which can cause poor approximation quality in finite element analysis. In ad- 
dition, meshes generated by isosurfacing techniques are often too dense. Therefore, im- 
proving mesh quality yet keeping the number of mesh elements small are two important 
issues that we will address in this section. 
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Figure 3 . Illustration of the surface generation and post-processing. (A) 
A 3D volume is first generated using the Gaussian blurring approach 
(equation (5.1)) from the molecule (PDB: 1CID). Shown here is part of 
the surface triangulation by the marching cube method. (B) The surface 
mesh after two iterations of mesh quality improvements. (C) After coars- 
ening, the mesh size becomes about seven times smaller than the original 
one. The mesh is also smoothed by the normal-based technique. 

5.2. Surface Mesh Improvement and Decimation. Surface mesh post-processing in- 
cludes quality improvement and mesh coarsening (decimation). The mesh quality can 
be improved by a combination of three major techniques: inserting or deleting vertices, 
swapping edges or faces, and moving the vertices without changing the mesh topol- 
ogy [19]. The last one is the main strategy we use to improve the mesh quality in our 
toolchain. For a surface mesh, however, moving the vertices may change the shape of the 
surface. Therefore, when we move the vertices, important features (e.g., sharp bound- 
aries, concavities, holes, etc.) on the original surface should be preserved as much as 
possible. To characterize the important features on the surface mesh, we compute so- 
called local structure tensor [18, 53, 57] as follows: 

/ M (*) (») (») (») (i) \ 
I rix' rix' n x riy rv x 'ny \ 



M' 



m = £ 



n , 



(5.2) 



\ (*) (*) W 0') (*) (0 I 

\ nl'nV nVn), ni'ny J 



<y 

where (n£ , riy \ n« ) is the normal vector of the i th neighbor of a vertex v and M' is 
the total number of neighbors. The normal vector of a vertex is defined by the weighted 
average of the normals of all its incident triangles. The local structure tensor basically 
captures the principal axes of a set of vectors in space. Let the eigenvalues of T(v) be 
Ai, A 2 , A 3 and Ai > A 2 > A 3 . Then the local structure tensor can capture the following 
features: (a) Spheres and saddles: Ai ~ A 2 ~ A 3 > 0; (b) Ridges and valleys: A x w 
A 2 > A 3 w 0; (c) Planes: Ai > A 2 w A 3 w 0. 

The quality of a mesh can be improved by maximizing the minimal angles. The angle- 
based method developed by Zhou and coauthors [61] utilizes this idea by moving a vertex 
(denoted by x) towards the bisectors of the angles formed by adjacent vertices on the 
surrounding polygon. This method works quite well for 2D planar meshes and has been 
extended in [55] for improving quadrilateral mesh quality as well. However, vertices 
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on a surface mesh can move with three degrees of freedom. If only the angle criterion 
is considered, the surface mesh may become bumpy and some molecular features may 
disappear. In other words, while the mesh quality is being improved, the geometric 
features on a surface mesh should be preserved as much as possible. To this end, we 
take advantage of the local structure tensor by mapping the new position x generated 
by the angle-based method to each of the eigenvectors of the tensor calculated at the 
original position x and scaling the mapped vectors with the corresponding eigenvalues. 
Let e 1; e 2 , e 3 denote the eigenvectors and Ai, A 2 , A 3 be the corresponding eigenvalues of 
the local structure tensor valued at x. The modified vertex x is calculated as follows: 



The use of eigenvalues as a weighted term in the above equation is essential to preserve 
the features (with high curvatures) and to keep the improved surface mesh as close as pos- 
sible to the original mesh by encouraging the vertices to move along the eigen-direction 
with small eigenvalues (or in other words, with low curvatures). Figure 3(B) shows the 
surface mesh after quality improvement, compared to the original mesh as shown in Fig- 
ure 3(A). Before quality improvement, the minimal and maximal angles are 0.02° and 
179.10° respectively. These angles become 14.11° and 135.65° after the improvement 
(two iterations). 

The surface meshes extracted by isocontouring techniques (e.g., the marching cube 
method) often contain a large number of elements and are nearly uniform everywhere. 
To reduce the computational cost, adaptive meshes are usually preferred where fine 
meshes only occur in regions of interest. The idea of mesh coarsening in our pipeline 
is straightforward — delete a node and its associated edges, and then re-triangulate the 
surrounding polygon. The local structure tensor is again used as a way to quantify the 
features. Let x denote the node being considered for deletion and the neighboring nodes 
be Vj, i = 1, • • • ,M, where M is the total number of the neighbors. The maximal length 
of the incident edges at x is denoted by L(x) = max^f 1 {c?(x, Vj)} where d(-, •) is the 
Euclidean distance. Apparently L(x) indicates the sparseness of the mesh at x. Let 
Ai(x), A 2 (x), A 3 (x) be the eigenvalues of the local structure tensor calculated at x, satis- 
fying Ai(x) > A 2 (x) > A 3 (x). Then the node x is deleted if and only if the following 
condition holds: 



where a and (3 are chosen to balance between the sparseness and the curvature of the 
mesh. In our experiments, they both are set as 1.0 by default. The threshold T is user- 
defined and also dependent on the values of a and (3. When a and (3 are fixed, larger T 
will cause more nodes to be deleted. For the example in Figure 3(C), the coarsened mesh 
consists of 8, 846 nodes and 17, 688 triangles, about seven times smaller than the mesh 
as shown in Figure 3(B). 

Mesh coarsening can greatly reduce the mesh size to a user- specified order. However, 
the nodes on the "holes" are often not co-planar; hence the re-triangulation of the "holes" 
often results in a bumpy surface mesh. The bumpiness can be reduced or removed by 
smoothing the surface meshes. We employ the idea of anisotropic vector diffusion [40, 
56] and apply it to the normal vectors of the surface mesh being considered. This normal- 
based approach turns out to preserve sharp features and prevent volume shrinkages [13] 




(5.3) 




(5.4) 
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better than the traditional vertex-based approach. Figure 3(C) shows the result after the 
mesh coarsening and normal-based mesh smoothing. 

5.3. Tetrahedral Mesh Generation. Once the surface triangulation is generated with 
good quality, Tetgen [44, 45] can produce tetrahedral meshes with user-controlled qual- 
ity. Besides the triangulated surface, our toolchain will have three other outputs for a 
given molecule: the interior tetrahedral mesh, the exterior tetrahedral mesh, and both 
meshes together. For the interior tetrahedral mesh, we force all atoms to be on the mesh 
nodes. The exterior tetrahedral mesh is generated between the surface mesh and a bound- 
ing sphere whose radius is set as about 40 times larger than the size of the molecule being 
considered. Figure 4 demonstrates an example of mesh generation on the mouse Acetyl- 
cholinesterase (mAChE) monomer. 




(A) (B) (C) (D) 

Figure 4. Illustration of biomolecular mesh generation. (A) The PDB 
structure of the mouse Acetylcholinesterase (mAChE) monomer. (B) The 
surface mesh generated by our approach. The active site is highlighted 
in yellow. (C) A closer look at the mesh near the active site. (D) The 
tetrahedral mesh between the molecular surface and the bounding sphere 
(not shown). 



6. Numerical Examples 

Two numerical examples with increasingly complexity of molecular surface are pre- 
sented to show the stability of the decomposition scheme and the convergence of the 
adaptive algorithm. In both examples, the Laplace equation for harmonic component is 
solved with finite element method. The gradient of the harmonic component is then com- 
puted and supplied for calculating the interface conditions of the regularized Poisson- 
Boltzmann equation. It is also possible to directly compute the harmonic component and 
its gradient from the solution representation for the Poisson equation of the harmonic 
component via surface integrals on the molecular surface. 

Example 1. The first numerical example is devoted to the comparison of two decom- 
position schemes discussed in the subsection 2. We use the model problem in Example 
2.1 because it admits an analytical solution for comparison. The computational domain 
is chosen to be a sphere with radius r = 5A. Figure 5 plots the computed regular poten- 
tial component and the full potential from the first decomposition scheme as well as their 
relative errors with respect to the analytical solutions, respectively. Chart B shows that 
the finite element solution of this regular component has an relative error below 3% over 
the entire domain. Because of the large magnitude of this regular potential, the absolute 
error is considerably large, see Chart A and in particular Chart C, where the analytical 
singular component is added to get the full potential. The amplification of the relative 
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error as analyzed in section 2 is seen from a comparison of Chart B and Chart D. This 
confirms that the first decomposition scheme is numerically unstable. 
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Figure 5. Solution of the Poisson-Boltzmann equation via the first de- 
composition scheme. (A): Computed regular component u r of the elec- 
trostatic potential (blue) versus the analytical solution (red). (B): Relative 
error in percentage of computed regular component u r of electrostatic 
potential. (C): Computed full electrostatic electrostatic potential (blue) 
versus the analytical solution (red). (D): Relative error in percentage of 
the computed full electrostatic potential. 

The numerical solutions via the second decomposition scheme demonstrate the desir- 
able numerical stability, as shown in Figure 6. The regular potential u r in Chart A is 
solved with the same mesh for Figure 5, and shows a very good agreement with the ana- 
lytical solution. The relative error is well below 1.5% over the entire domain, and is well 
below 0. 1% in the interior of the molecule. Compared to Figure 5, it is seen that the mag- 
nitude of the regular component of the stable decomposition is much smaller. Because 
the harmonic and regular components are both solved numerically, it is worthwhile to 
examine the summation of these two numerical solutions and compare the total with the 
exact solution; this is plotted in Chat B. The discontinuity indicates that the decompo- 
sition is only applied inside the biomolecule, and that the harmonic component is much 
larger than the regular component. This further suggests that the overall relative numeri- 
cal error inside the biomolecule will be larger than that in the solvent region. Interesting 
enough, most intermolecular electrostatic interactions are occurred through the solvent, 
and thus the stable decomposition can still provide the electrostatic potential of high fi- 
delity for describing these interactions. We then refine this mesh globally by bisecting all 
the edges; the relative error is reduced 0.5% in most of the domain except in the vicinity 
of the dielectric interface where the error does not show a noticeable decrease, see Chart 
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D. This is because the middle point of an edge on the interface maybe not located on 
the interface and therefore violates the assumption on the discretization of the molecular 
surface, thus the approximations to the interface and to the interface conditions are not 
improved with this globally refinement. To satisfy this assumption we apply this global 
refinement first and then move the middles points of all the interface edges back to the 
interface. This new refinement approach successfully scales down the numerical error 
near the interface, see Chart E. 

Both two decomposition schemes give rise to an elliptic interface problem, whose 
solution is of C° only and can be appreciated from Chart A of Figures 5 and 6. 

Example 2. This second numerical experiment is conducted on an insulin protein 
[51] (PDB ID: 1RWE). This protein has two polypeptide chains, one has 21 amino acid 
residues and the other has 30 residues, and has 1578 atoms in total. Because there is no 
analytical solution available for accuracy assessment we solve the Poisson-Boltzmann 
equation on four progressively refined meshes and use the solution on the finest mesh as 
the reference to measure the accuracy of other three solutions. In Table 1 we show the 
computed electrostatic solvation energy (AG eie ) and the corresponding relative error in 
the solution (eAc? eie ) f° r eacn solution. This energy is defined as 

AC7 e 2 e = - / (psoi - p vac )p f dx, 

where p so i is the electrostatic potential of the solvated molecule while p vac is the potential 
for the molecule in vacuum, where the dielectric constant is assumed to the same as the 
interior of the molecule. It turns out that p vac is essentially the singular component inside 
molecule, and thus the solvation energy can be directly computed as 

AG eZe = I [ (p h + p r )p f dx. 

Assume that the solution on the finest mesh is convergent, we computed the relation 
error in the solution energy for three coarser meshes. The diminishing of this relative 
error confirms the convergence of the numerical method for computing electrostatics of 
realistic biomolecules. 



level of mesh 


# of tetrahedra 


# of nodes 


# of nodes on T 


AG e i e 


e AG eie 


1 


280928 


45119 


7681 


-1476.9 


0.0757 


2 


340410 


54685 


10257 


-1403.6 


0.0224 


3 


474011 


76036 


14982 


-1380.3 


0.0054 


4 


626221 


100393 


19816 


-1372.9 





Table 1 . The electrostatic solvation energy AG e i e and the correspond- 
ing relative error e& Gelc for progressively refined meshes, confirming con- 
vergence of the discretization technique based on the new regularization. 



7. Summary 

In this article, we considered the design of an effective and reliable adaptive finite 
element method (AFEM) for the nonlinear Poisson-Boltzmann equation (PBE). In Sec- 
tion 2, we began with a very brief derivation of the standard form of the Poisson- 
Boltzmann equation. We examined the two-scale regularization technique described 
in [11], and briefly reviewed the solution theory (a priori estimates and other basic re- 
sults) developed in [1 1] based on this regularization. We then described a second distinct 
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Figure 6. Solution of the Poisson-Boltzmann equation via the second 
decomposition scheme. (A): Computed regular component u r of the elec- 
trostatic potential (blue) versus the analytical solution (red). (B): Com- 
puted regular component u r plus the harmonic component u h of the elec- 
trostatic potential (blue) versus the analytical solution (red). (C): Relative 
error in percentage of the computed regular component of the electrostatic 
potential on an initial mesh. (D): Relative error in percentage of the com- 
puted regular component of the electrostatic potential; globally refined 
mesh. (E): Relative error in percentage of the computed regular com- 
ponent of the electrostatic potential; mesh locally refined on molecular 
surface and the boundary. 

regularization and explained why it is superior to the original approach as a framework 
for developing numerical methods. We then quickly assembled the cast of basic mathe- 
matical results needed for the second regularization. In Section 3, we described in detail 
an adaptive finite element method based on residual-type a posteriori estimates, and 
summarized some basic results we needed later for the development of a corresponding 
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Figure 7. Left: The electrostatic potential mapped on the molecular 
surface of the insulin protein. Right: The surface mesh of insulin protein 
in the finite element model. 



convergence theory. We presented this new convergence analysis in Section 4, giving the 
first AFEM contraction-type result for a class of semilinear problems that includes the 
Poisson-Boltzmann equation. 

We gave a detailed discussion of our mesh generation toolchain in Section 5, includ- 
ing algorithms designed specifically for Poisson-Boltzmann applications. These algo- 
rithms produce a high-quality, high-resolution geometric model (surface and volume 
meshes) satisfying the assumptions needed for our AFEM algorithm. These algorithms 
are feature-preserving and adaptive, designed specifically for constructing meshes of 
biomolecular structures, based on the intrinsic local structure tensor of the molecular 
surface. Numerical experiments were given in Section 6; all of the AFEM and mesh- 
ing algorithms described in the article were implemented in the Finite Element Toolkit 
(FETK), developed and maintained at UCSD. The stability advantages of the new regu- 
larization scheme were demonstrated with FETK through comparisons with the original 
regularization approach for a model problem. Convergence and accuracy of the AFEM 
algorithm was also illustrated numerically by approximating the solvation energy for a 
protein, in agreement with theoretical results established earlier in the paper. 

In this article, we have examined an alternative regularization which must be used in 
place of the original regularization proposed in [11], due to an inherent instability built 
into the original regularization. We showed that an analogous solution and approxima- 
tion theory framework can be put into in place for the new regularization, providing 
a firm foundation for the development of a large class of numerical methods for the 
Poisson-Boltzmann equation, including methods based on finite difference, finite vol- 
ume, spectral, wavelet and finite element methods. Each of these methods can be shown 
to be convergent for the regularized problem, since it was shown in this article to allow 
for a standard H 1 weak formulation with standard solution and test spaces. Our primary 
focus in this article then became the development of an AFEM scheme for the new reg- 
ularized problem, based on residual-type a posteriori error indicators, a fairly standard 
and easy to implement marking strategy (Dorfler marking), and well-understood simplex 
bisection algorithms. We showed that the resulting AFEM scheme is reliable, by proving 
a contraction result for the error, which established convergence of the AFEM algorithm 
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to the solution of the continuous problem. The AFEM contraction result, which is one 
of the first results of this type for nonlinear elliptic problems, follows from the global 
upper boundedness of the estimator, its reduction, and from a quasi-orthogonality re- 
sult that relies on the a priori L°° estimates we derived. This new AFEM convergence 
framework is distinct from the analysis of nonlinear PBE with the previous regulariza- 
tion approach from [1 1], is more general, and can be applied to other semi-linear elliptic 
equations [27]. The contraction result creates the possibility of establishing optimality 
of the AFEM algorithm in both computational and storage complexity. 

We note that for computational chemists and physicists who rely on numerical solu- 
tion of the Poisson-Boltzmann equation, discretizations based on the stable splitting as 
described in the current paper are the only reliable numerical techniques under mesh 
refinement for the Poisson-Boltzmann equation that we are aware of (both provably con- 
vergent and stable to roundoff error). While one must take care with evaluation of the 
singular function u s , since this generally involves pairwise interactions between charges 
and mesh points, the alternative to using these types of splitting discretizations is to lose 
reliability in the quality of the numerical solution. While we focused on (adaptive) finite 
element methods in this article, we emphasize that the splitting framework can be eas- 
ily incorporated into one's favored (finite difference, finite volume, spectral, wavelet, or 
finite element) numerical method that is currently being employed for the PBE. 
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